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Abstract In this paper, we examine linear programming (LP) relaxations 
based on Bernstein polynomials for polynomial optimization problems (POPs). 
We present a progression of increasingly more precise LP relaxations based on 
expressing the given polynomial in its Bernstein form, as a linear combination 
of Bernstein polynomials. The well-known bounds on Bernstein polynomials 
over the unit box combined with linear inter-relationships between Bernstein 
polynomials help us formulate “Bernstein inequalities” which yield tighter 
lower bounds for POPs in bounded rectangular domains. The results can be 
easily extended to optimization over polyhedral and semi-algebraic domains. 
We also examine techniques to increase the precision of these relaxations by 
considering higher degree relaxations, and a branch-and-cut scheme. 

Keywords Polynomial Optimization Problem • Bernstein Polynomials • 
Linear programming 


1 Introduction 

In this paper, we examine linear programming relaxations for polynomial op¬ 
timization problems (POP) that seek to optimize a multivariate polynomial 
p(x) over a compact interval domain x S [£, m]. Our approach is based on two 
ideas: (a) We consider a reformulation of the problem as a linear program 
using Bernstein polynomials. However, doing so also increases the number 
of decision variables and constraints in the problem, (b) Next, we present 
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valid inequalities for improving the approximation. These inequalities are de¬ 
rived from well-known properties of Bernstein polynomials that yield linear 
inter-relationships between the decision variables of the linear program. Our 
approach is extended to handle compact domains described by semi-algebraic 
constraints. We present a branch-and-cut scheme that introduces the cutting 
plane inequalities hand-in-hand with a decomposition of the feasible region. 

The problem of optimizing polynomials over an interval is well-known to 
be non-convex, and is in fact NP-hard. Nevertheless, well-known classes such 
as linear, quadratic, or even integer linear programs can be viewed as partic¬ 
ular cases of POPs. Also, since polynomials provide a good approximation for 
non linear functions, solving POPs efficiently is a big step toward handling 
more complex problems. Finally, a lot of problems arising from disparate do¬ 
mains such as biology, robotics and engineering can be formulated as POPs. 
Our interest is motivated by verification and synthesis problems for dynamical 
systems such as safety, reachability and stability verification. These problems 
can be reduced to POPs. In fact, the motivation of this paper comes from 
our previous work, where we aim to prove stability for polynomial dynamical 
systems m- Therein, Bernstein polynomials were used as an alternative to 
the well-known sum of squares (SOS) approach in order to avoid the numer¬ 
ical issues of semi-definite programming (SDP) [T^ITlfTB]. In this regard, the 
Simplex algorithm can be implemented in exact arithmetic to yield numeri¬ 
cally validated lower bounds to the optimal value of the POP, thus formally 
establishing the stability of the process. The success of the approach in a large 
set of benchmarks motivates us to go further, improve the results and make 
them known in an optimization context. 

More precisely, we show how POPs can be relaxed to linear programs 
thanks to the use of Bernstein polynomials, and a well-known reformulation- 
linearization technique (RLT) described by Sherali et al [TBIII7] . In fact, the 
properties of Bernstein polynomials inside the unit box offer us an elegant ap¬ 
proach to improving the RLT approach. We formulate these properties as linear 
inequalities to obtain guaranteed lower (upper) bounds for our minimization 
(maximization) problems. This will be useful in cases where the POP does not 
need to be solved exactly. In the latter case, we combine our inequalities with 
a branch-and-bound decomposition process originally described by Nataraj et 
al [ig. 

We evaluate our approach using a set of benchmarks described in Nataraj 
et al m to characterize the effect of adding the extra Bernstein inequalities 
to the RLT approach. We observe that while the addition of these inequali¬ 
ties improves the lower bound, it is not sufficient for yielding tight bounds. 
Next, we consider the addition of Bernstein inequalities in a “branch-and-cut” 
approach that combines the addition of cutting planes “on-demand” with a 
branch-and-bound decomposition of the domain. We find that all approaches 
eventually yield tight bounds on the value of the global optimum. Therefore, 
we compare the computational time for various approaches. Finally, we com¬ 
pare the various approaches on benchmarks from our previous work m which 
consists of a set of polynomial Lyapunov functions used as stability proofs for 
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polynomial dynamical systems. In this particular case, our goal is to show that 
the functions are non-negative over a domain. We adapt the branch-and-bound 
scheme for this application to evaluate its effectiveness. 

The results of our evaluation are mixed: we observe that adding cutting 
plane inequalities does result in tighter lower bounds on the optimum and 
therefore examining fewer cells in the branch-and-bound approach. However, 
this comes at the cost of obtaining larger linear programs due to the extra 
inequalities, and therefore, an overall larger computation time. We show that 
the careful consideration of inequalities to be introduced yields a “sweet spot” 
for better approximations using less computation time. 


1.1 Organization 

In Section [21 we present basic notions and properties related to Bernstein 
polynomials. All the results of this section are quite standard , therefore proofs 
are omitted. Section [3] is the core of the paper. In this section, Bernstein 
polynomials and their properties inside the unit box are translated into a 
series of inequalities yielding a corresponding set of LP relaxations of increasing 
precision. An iterative approach mixing these relaxations is presented, and a 
criterion for checking if the given lower bound meets the optimal value of the 
original problem are also given. In Section 01 we show how bounds can be 
made arbitrary tighter using some techniques such as decomposition (branch- 
and-bound scheme). 


1.2 Related Work 

Since solving a POP is generally NP-hard, existing work consists of relaxing 
it in order to obtain an easier problem for which efficient solvers exist. In the 
literature, we can distinguish two types of relaxations. The first class is called 
LP relaxations. These approaches approximate the POP using linear programs 
that can be efficiently solved using an LP solver. A popular LP relaxation is the 
reformulation linearization technique (RLT) given by Sherali et al [16U17] . The 
approach was improved by Nataraj et al m for solving POPs, wherein the use 
of the Bernstein basis was proposed as an improvement. In particular, Nataraj 
et al made use of the property that Bernstein polynomial coefficients over a box 
form a lower bound of the polynomial. In this work, we show that this property 
is simply the optimal value of a LP formed by a series of inequalities that relate 
one Bernstein polynomial to another. In doing so, we formulate numerous valid 
inequalities that improve substantially on this bound. Another recent approach 
called DSOS (Diagonally-dominant Sum of Squares) was formulated by Ali 
Ahmadi et al [T] by relaxing positive semi-definiteness of a matrix using the 
stronger condition of diagonal dominance. In fact, Ali Ahmadi’s approach can 
be seen as selecting a finite set of generators from the infinitely generated cone 
of positive polynomials in the polynomial ring R[x]. In contrast, our approach 
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also adds a finite set of generators to the cone of positive polynomials over 
a compact interval. Naturally, both choices of finite bases involve a tradeoff 
that are optimal for certain classes of problems. In particular, we choose the 
Bernstein polynomials and utilize the set of linear inter-relationships between 
these. Extending our approach to possibly cover the polynomial basis used in 
the DSOS approach is currently under investigation. 

As an alternative to LP relaxations, we can formulate SDP relaxations. 
In 2001, Lassere proposed what was called a Linear matrix equality (LMI) 
relaxation [7]. The main idea is to map the polynomial optimization problem 
to an optimization problem over probability measures and then use results 
from moment theory. Subsequently, Parillo introduced the SOS programming 
approach that has become one of the most popular SDP relaxations [13] . The¬ 
oretically, following the comparison made by Lasserre [5] between SDP (LMI) 
and LP (RLT) relaxations, one concludes that the SDP approach is much more 
precise at the extra (polynomial) cost of solving an SDP. In fact the compar¬ 
ison points out that for the LP (RLT) relaxation, convergence results to the 
optimal value are not always guaranteed, in contrast to SDP relaxations. Also, 
the comparison shows that RLT cannot be exact whenever the global optimum 
belongs to the interior of the feasible set. We will show in this paper, that this 
claim does not remain true (see Example [2]) when the Bernstein inequalities 
suggested here are used. Furthermore, in practice, the SDP approach suffers 
from numerical issues. This was pointed in our previous work m when us¬ 
ing SOS programming for Lyapunov function synthesis. Other approaches like 
interval methods uni and decomposition techniques exists. In this paper, we 
will focus on the related scheme given by Nataraj et al in m, since it is fully 
based on the use of Bernstein coefficients. We will build on this approach by 
adding the extra Bernstein inequalities. 


2 Overview of Bernstein Polynomials 

Bernstein polynomials were first proposed by Bernstein as a constructive proof 
of Weierstrass approximation theorem |3] , and are useful in many engineering 
design applications for approximating geometric shapes |5] . They form a basis 
for approximating polynomials over a compact interval, and have nice proper¬ 
ties inside the unit box (see [TT] for more details). We first examine Bernstein 
polynomials and their properties for the univariate case, and then extend them 
to multivariate polynomials (see [HIS])- 

Definition 1 (Univariate Bernstein Polynomials) Given an index 
i € {0,..., m}, the univariate Bernstein polynomial of degree m over [0,1] 
is given by the following expression: 



( 2 . 1 ) 
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Using these polynomials, monomials can be written as follows: 


X = 


j^i 


m \ 
*) 


for alH = 0,..., m. 


( 2 . 2 ) 


Then, in the Bernstein polynomial basis, polynomial p{x) 
degree m can be written as: 


m 


of 

i=o 


m 

z=0 


where for alH = 0,..., m: 




(2.3) 


The coefficients bi^rn are called the Bernstein coefficients of the polynomial p. 

Bernstein polynomials have many interesting properties on the unit interval 
[0,1]. We summarize the most relevant ones for our applications. 

Lemma 1 Bernstein polynomials have the following properties: 

m 

1. Unit partition: = 1. 

2—0 

2. Bounds: 0 < ffimix) < Vi = 0,... ,m. 

3. Induction: Vi = 0,..., m - 1. 

Using these properties, the following result holds: 

Corollary 1 On the interval [0,1], the following inequality holds 


min hi^rn<p{x)< max b^^rn- (2.4) 

2=0,...,m 2=0,...,m 


The equality min bim= min p{x), respectively max bim= max p{x), 

’ a::€[0,l] 2 = 0 ,...,m ’ a;G[0,l] 

holds iff min bi^rn & {bo.m,bm,m}, respectively max G {6o,m, 

2=0,...,m 2=0,...,m 

This is commonly called the vertex condition. 

We generalize the previous notions to the case of multivariate polynomials 
i.e p{x) = p{xi,... ,Xn) where x = (xi,..., Xn) € U = [0, Ij”. For multi 
indices, I = (ii,... ,i„) G N", J = {ji,... ,j„) G N”, we will use the following 
notation throughout this paper: 

“ I + J = {il + jl, ■ ■ ■ ,in + jn)- 
— X^ = xf^ X X2^^ • • • X Xn^’'. 
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— I < J ii < ji, for alH = 1,... n 



— Ir,k = (*i, • ■ •, + k, ir+1, ■ • ■ , in) where re {1,..., n} and k Gh. 

Let us fix our maximal degree <5 = (di,..., 5n) G N” for a multivariate polyno¬ 
mial p {Si is the maximal degree of xi for alH = 1,..., n). Then the multivariate 
polynomial p can be written as: 

p{x) = E pix^ where pi G K, V/ < S. 

I<5 


Multivariate Bernstein polynomials are given by products of the univariate 
polynomials: 


Bps{x) = Pi^,s^{xl)...|3^„,s„ixn) = [ / ]xj'{l-Xj)^^ 

(2.5) 

Thanks to the previous notations, these polynomials can also be written as: 


5-1 


Bps{x) = (^/j ^ (1" “ 

The expression of monomials using these polynomials is: 


( 2 . 6 ) 


= E 


I<J<5 


Bj^s{x), for all / < (5 


(2.7) 


Now, we can give the general expression of a multivariate polynomial in 
the Bernstein basis: 

Pix) = y^,bi^sBi^s{x), 

I<5 

where Bernstein coefficients {bi^s)i<s are given as follows: 


j<i 



( 2 . 8 ) 


Therefore, the generalization of Lemma [T] will lead to the following properties: 


Lemma 2 For all x = (cci,..., Xn) G U we have the following properties: 


1. Unit partition: E Bpsix) = 1. 
i<s 


2. Bounds: 0 < Bj s{x) < Bj^s{^), for all I < 6. 

3. Induction: BpSr,-i = Vr = 1, ..., n. 


I < Sr,-1. 
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Also, Corollary [T] can be generalized as follows: 

Corollary 2 Let p be a multivariate polynomial of degree 5 over the unit box 
U = [0, 1]" with Bernstein coefficients bj^s where I < S. Then, for all tc G U, 
the following inequality holds: 

mini)/ s < p(x) < max6/ s- (2.9) 

i<6 ’ i<s 

The vertex condition holds iff the minimum value (respectively the max¬ 
imum value) is reached for an index I* G Sg where: 

Sg = {I = (ii, ... ,in) G N", such that ij G {0, Sj}, Vj = 1,...,n}. 

Given the Bernstein coefficients (&/, 5)/<5 for a polynomial p, the vertex 
condition is quite easy to check using the steps outlined below: 

1. Find I* := argm.mj^g{bi ^s)■ 

2. Check for each j e [i, n] if I* = 0 or I* = 6j. 

3. If the previous step succeeds, vertex condition holds and bj^^s is a global 
minimum of p inside the unit box. Otherwise, vertex condition fails. 

Checking the vertex condition will be an important primitive for the overall 
approach that will be developed in this paper. 

Finally, consider an arbitrary, bounded interval 1C : [£]_,^] x • • • x \xn,^], 
wherein —oo < Xj < xj < oo, for all j = 1,..., n. It suffices to a map K, into 
the unit box U by applying the following change of variables from x to z: 
Zj = — for all j = 1,..., n. Doing so, the results from Lemma [2] can be 
transferred to arbitrary boxes 1C. 


3 Bernstein Polynomial Relaxations for Polynomial Optimization 
Problems 

Given a multivariate polynomial p and a rectangle 1C, we consider the following 
optimization problem : 

minimize p{x) , , 

s.t xGK.. 

Whereas dm is hard to solve, we will construct a linear programming (LP) 
relaxation, whose optimal value is guaranteed to be a lower bound on p*. In 
this section, we will use Bernstein polynomials for the unit box (/C = [0,1]"). 
If /C is a general rectangle, we use an affine transformation to transform p and 
)C back to the unit box. 




Mohamed Amin Ben Sassi, Sriram Sankaranarayanan 


3.1 Reformulation Linearization Technique (RLT) 

We first recall a simple approach to relaxing polynomial optimization problems 
to linear programs, originally proposed by Sherali et al. [MHZ]. We then carry 
out these relaxations for Bernstein polynomials, and show how the properties 
in Lemma [5] can be incorporated into the relaxation schemes. Recall, once 
again, the optimization problem (13.11) over the unit box /C: 

p* = minimize p{x) 

s.t a; G/C =: [0,1]", 

n 

where fC is represented by the constraints K : /\ Xj > 0 A {1 — Xj) > 0. The 

i=i 

standard RLT approach consists of writing p{x) = 'Y^jPix^ as a linear form 
p{x) : 'YhiPiyi for fresh variables yi that are place holders for the monomials 
x^. Next, we write down as many facts about x^ over /C as possible. The basic 
approach now considers all possible power products up to a maximal degree D 
i.e of the form ttj^s '■ x'^ {\—x)^~'^ for all J < 5 where |5| = D. Clearly if a; G /C 
then TTj^six) > 0. Expanding ttj^s in the monomial basis as ttj^s ■ J2i<s o-i,jx^^ 
we write the linear inequality constraint 

ai,,j yi > 0. 

i<j 

The overall LP relaxation is obtained as 
minimize E piyi 

E (3.2) 

s.t. 2, ^i.J yi ^ Oj for each J < 6. 

i<j 

Additionally, it is possible to augment this LP by adding inequalities of the 
form ij < yi < uj through the interval evaluation of over the set 1C. 

Remark 1 The extra “facts” that form the constraints in Eq. 13.21 are akin to 
valid inequalities or cuts that incrementally refine an over-approximation of the 
feasible region. Unfortunately, the number of such inequalities is exponential in 
|(5|. Rather than adding these all at once to yield a single LP, we may add them 
on demand, iteratively solving a series of LPs wherein the new inequalities are 
introduced as cutting planes to help improve the solution. 

Proposition 1 For any polynomial p, the optimal value computed by the 
LP ()3.2|) is a lower bound to that of the polynomial program en). 

Example 1 We wish to solve the following POP (or find a lower bound for its 
solution): 


minimize Xi^ + xf^ 
s.t. (a;i,a; 2 ) G [0,1]^ 


(3.3) 
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Using the RLT technique for a degree D = 2 we denote by yij the fresh 
variables replacing the non linear terms = a;i®a; 2 -^ for all {i,j) G such 
that i + j <2. We obtain an LP which is shown, in part, below: 

minimize 3/2,0 + 2/0,2 
s.t 

0 < 2/2,0 < 1 
0 < 2/0,2 < 1 


The optimal solution obtained from the LP is 0, which coincides with the 
optimum of the original problem. 


3.2 RLT using Bernstein Polynomials 

The success of the RLT approach depends heavily on writing “facts” involving 
the variables y/ that substitute for We now present the core idea of using 
Bernstein polynomial expansions and the richer bounds that are known for 
these polynomials from Lemma [2] to improve upon the basic RLT approach. 

Linear rel 2 Lxations : First, we write p(x) as a weighted sum of Bernstein 
polynomials of degree (5. 

p(x) : bpsBps , 

i<s 

wherein bj s are calculated using the formula in equation (12.81) . Let us introduce 
a fresh variable zj^s as a place holder for Bi^s{x). Lemma d] now gives us a 
set of linear inequalities that hold between these variables zps- We formulate 
three LP relaxations, each providing a better approximation for the feasible 
region of the original problem (EH). 

= minimize bj^gzj^s 
i<s 

s.t. ^,>0 I <5, 

= 1 

I<S 

zi^s s R, I < S. 

Remark 2 It is easy to see that pf''^ = min/i /_5 (the smallest Bernstein coef¬ 
ficient). As a result, it can be computed quite efficiently without actually in¬ 
voking an LP solver. In fact, the branch-and-bound approach of Nataraj m 
is based on this relaxation. 
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Using the upper bound on Bernstein polynomials from Lemma [21 we can 
strengthen (13.411 further, as follows: 

= minimize 6/,5-Z/,5 
i<s 

s.t >0, I < S, 

zi,s < Bj^s (I) I < S, -(r- Upper Bounds (3.5) 

1 

I<S 

zij G R, I < 6. 

Next, tighter relaxation can be obtained by adding the induction relations 
between Bernstein polynomials of lower degrees. More precisely, using in ad¬ 
dition the third property of Lemma |21 we obtain the following linear program: 

= minimize bi^szps 
I<5 

s.t zi,K G M, I < K, K < S, 

0 < zi.K < I<K, K<6 

zi,K = 1 , K <s 

I<K 

zi.Kr--! = z:i,k + Vr e {1,... ,n} s.t /^.i < K. 

(3.6) 

Remark 3 Note that Eq. (13.5p involved variables zj^s for I < 6. The formu¬ 
lation in Eq. (|3ll) involves a larger set of “lower degree” terms of the form 
zi,K wherein I < K and K < 5. These terms are, in fact, not necessary as 
demonstrated in Prop. [31 

Each of these relaxations provides a lower bound on the original polynomial 
optimization problem. 

Proposition 2 < ps^^^ < Ps^'^^ < P* ■ 

Proof We already know thanks to Corollary |T] that ps^^'* < p*. 

Now, consider any feasible solution y to the problem 6ID which is equiv¬ 
alent to 

minimize E bpsBj^siy) 
i<s 

s.t ye [0,1]". 

We note that replacing zj^g = Bj^siv) the vector of all zi^g form a feasible 
solution to each of the two relaxations. Therefore, p^p < p* for all j S {1,2}. 

Also it is easy to see that these relaxations are increasing (since they are 
constructed by adding extra constraints). Therefore, ps^°^ < ps^^'^ < Ps^"^^ ■ 

□ 
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Example 2 Let’s consider p{x) = on [—1,1]. For a degree 5 = 2, the min¬ 
imum of Bernstein coefficient is = —1. Whereas using (E3|), we found 
ps^^^ = 0 which coincides with the optimum. 

Now, we consider the polynomial p{x,y) = x'^ + on [—1,1]^, we found 
ps^^^ = —2 and ps^^'' = —0.5. Using (13.61) . we obtain ps^"^^ = 0 which is the 
exact optimal value. 

Now, in order to simplify relaxation dSH), we formulate an equivalent relax¬ 
ation that only uses decision variables This is achieved by replacing lower 
degree variables {zpK where K < 5) hy a matrix product involving variables 
zps ■ More precisely, we have the following result : 

Proposition 3 There exist a matrix As and a veetor cs such that the LP 
formulation in Eg. (EH) can be written as 

ps^^'^ = minimize bs ■ zs 

S.t O 5 ^ Zs ^ Us, /Q 

Is ■ zs = 1, ^ ^ 

Aszs < cs, 

wherein the notation as stands for a vector {ai^s)i< 5 , uj^s = Bj^siy); 0i,s = 0 
and = 1. 

Proof Each Bernstein polynomial Bpxix) can be written uniquely as 

Bpk{x) = ^ b[j'^^Bj^s{x) 
j<s 

wherein J < 5 form the Bernstein coefficients for the polynomial 

Bi^k{x). Translating this, we obtain the relation 

zi,K = X] b^j'^'^zys = ■ zs 

J<5 

wherein is the vector of Bernstein coefficients (b^j’^^)j<6 and zs = 

{zj,s)j< 5 - The result can now be established by systematically replacing each 
variable zj^k for K < 5 in (13.611 into an expression in terms oi zs- □ 

Remark f The computation of the pair {As,cs) in Eq. (13.711 depends only on 
<5, and is independent of the actual objective function. As a result, it can 
be computed offline, once for a given problem setup in terms of number of 
variables and <5. 

Iterative approach: In many cases, the optimal value given by the linear 
program EH can be obtained with fewer number of constraints i.e instead of 
having the constraints given by the pair {As, cs) only some of them are needed. 
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2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


Data: POP objective p{x), Constraints gi{x) < 0,... ^gKix) < 0, Box B, Degree 
limit vector 5. 

Result: p|: an underapproximation to the POP. 
begin 

Transform p, 51 ,... , over the unit box K : [0,1]^ 

Compute matrices {As,cs) 

Initialize (As^cs) empty matrices 
changeOccurred := TRUE 
while changeOccurred do 

:= Solution to the LP (13.811 
changeOccurred := FALSE 
for each row j in (A^, cs) do 
if Asjzs > Cs then 

changeOccurred := TRUE 
Add row j from Aj, eg to Ag , eg 
Remove row j from {Ag, cg) 


Algorithm 1: Overall algorithm for solving a POP using iterative Bernstein 
polynomial relaxation. 


In fact, often a large number of constraints are inactive for the optimal solution. 
More precisely, we solve LPs of the form; 


= minimize bg • zg 

S.t Oi < 05 < ug, 

Ig ■ zg = 1 , 
Agzg < Cg, 


(3.8) 


where {Ag, eg) contains a subset of the rows in the matrix {Ag, eg). Algorithm[T] 
shows the overall iterative scheme. 

1. Lines[2]to |4]show the initialization steps that involve computing the matri¬ 
ces {Ag,cg). The incremental computation involves using matrices (A 5 ,ci) 
that are initially empty. 

2. Solve the linear program p.8l) . which is initially the same as (13.51) . 

3. At each step, we obtain the current optimal value pg and an optimal solu¬ 
tion Zg. 

4. The for loop in line [13] iterates through all rows j of the matrix Ag such 
that Agjzg < cgj is violated. 

5. We these violated rows to the linear program (13.5L remove them from 
{Ag, eg). 

6. Termination happens whenever no violated rows are found in the for loop. 

Exact relaLxation: The decision variables zj^g introduced during the RLT 
technique are fresh variables that substitute the nonlinear polynomials Bi^g{x). 
A sufficient condition for an exact relaxation to hold is that optimal solutions 
zi^g* = Bpg{x*) for all / < (5 where x* € 1C. It is easy to see that when this 
happens, x* is in fact a global optimum for our problem. 
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Proposition 4 Let be the optimal value given by our relaxation. If there 
exist X* S /C(:= [0,1]") such that 

(A 

= E (3.9) 

0<J<<5 \l) 

then the relaxation is exact i.e p* = ps* and x* is the global minimum. 

Proof The conditions (13.9p can be written as a BgZg = {x*^)i<s where Bs is 

the matrix given by (13.91) . If this condition holds then we have 

for all / <6 which implies that ps* = ’'^^bi^sBj^six*) = p{x*). This shows 

I<5 

that Ps* is also an upper bound and prove that the condition is sufficient. 

The converse is not necessarily true: it is easy to construct examples 
wherein is optimal andp| coincides with a global optimum, but Zg ^ Bsfx*) 
for any x* in the domain. In fact, since Zg is not unique, then we can have 

zi,s* Bps{.x*) whereas '^bj^szps* = '^bpsBps{x*). 

i<s i<s 

Given an optimal solution Zg, we now provide a procedure that attempts 
to possibly find a x* G 1C such that Bsfx.*) = Zgi 

1. Each variable Xj is itself a polynomial and thus can be written uniquely in 

the Bernstein form as Xj : J2i<s Bi^s{x), wherein aP are the Bernstein 
coefficients of Xj . 

2. Therefore, compute a nominal vector x as Xj : ^jcg aPzf g. 

3. Use X to check if Zg = Bs(x). If yes, we conclude exactness of our relaxation 
with X* = X, and stop. Otherwise, we conclude that no such x* exists. 

Example 3 Consider the problem in Example [21 For the univariate case, one 
can check that relaxation (Ell) is exact. In fact, the optimal solution is = 
(0.25,0.5,0.25). Using the previous remark, we found x = 0.5 and we check 
that zps* = Bps{x) for all I < S. Then the relaxation is exact and x* = 0.5 
which corresponds to zero after a linear transformation to [—1,1]. For the 
bivariate case, relaxation (13.61) is exact but the condition (13.91) does not hold. 
This is due to the fact that is not unique. 


3.3 Numerical examples 

Thus far, we have presented three LP relaxations using Bernstein polynomials. 
For the formulation in Eq. (ESI), we provide a technique to reduce the number 
of variables by computing matrices {As,cs) that substitute constraints over 
variables zjioi K < <5 in terms of variables zs (Eq. (13.71) 1. Next, we provide 
an iterative approach that avoids an upfront solution to Eq. (1X71) . considering 
an iterative and incremental addition of constraints as in Eq. (IXX) . Also, our 
approach thus far is monolithic: we translate a single instance of a POP into 
a LP formulation without considering subdivisions of the feasible region 1C. 
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Table 1 Performance of the “monolithic” Bernstein relaxation for benchmark problems 
taken from Nataraj et al m Legend: ID: problem ID as given in [S], 5: the maximum 
degrees on variables, Minimum Bernstein coefficient, LP relaxation in Eq. 113.5l l. 

p^^^: LP relaxation in Eq. 113.711 . #r(A 5 ): number of constraints for formulation in Eq. 113.711 . 
^y{As): number of rows for the reduced problem in Eq. 113.811 . k: number of iterations, t\ 
time (seconds) to compute the matrices 


ID 

5 





#r(Aj ) 

k 

t 

1 

(4,4) 

-1170 

-911.47 

-856.42 

200 

6 

3 

0.1 

2 

(6,4) 

-7990.8 

-7195 

-6709.9 

385 

3 

2 

0.2 

3 

(2,2) 

-926 

-451 

-316 

27 

4 

2 

0.1 

4 

(4,2) 

-9994 

-6223.4 

-4721.4 

75 

6 

1 

0.1 

5 

(2,2,2) 

-240 

-109.5 

-66.75 

189 

9 

1 

0.1 

6 

(2,4,4) 

-200299 

-199930 

-139355.28 

1563 

28 

2 

1.4 

7 

(1,2,1) 

-36.7127 

-36.7127 

-36.7127 

42 

0 

1 

0.1 

8 

(2,4,4) 

-20218 

-19948 

-14290.38 

1692 

35 

3 

1.4 

9 

(1,1,3,3) 

-3.77 

-3.77 

-3.53 

836 

4 

1 

0.8 

10 

(1,2,2,2) 

-25.2 

-21.35 

-21.35 

594 

0 

0 

0.5 

11 

(2,2,2,2) 

-1020 

-542 

-260 

2700 

36 

1 

1.7 

12 

(1,1,1,1,2) 

-55 

-50 

-32.5 

438 

8 

2 

0.4 

13 

(2 ,...,2) 

-11 

-6.58 

-0.5 

45927 

449 

2 

1630 

14 

(1,2,2,3,1,1) 

-1.44 

-1.44 

-1.44 

9432 

0 

0 

64.8 

15 

(2 ,...,2) 

-13 

-7.5 

- 

- 

- 

- 

TO 

16 

(2 ,...,2) 

-2.04 

-2.02 

- 

- 

- 

- 

TO 


We evaluate these techniques using benchmark examples proposed by Nataraj 
et ah m- In TablelU we report the optimal values of the proposed relaxations, 
the size of matrices As ^ As, the number of iteration needed and the compu¬ 
tation time for the matrix As (we print ‘TO’ if the computation time exceeds 
30 minutes). We find that considerable reduction is made by considering As 
instead of As and also a considerable improvement in the lower bound is ob¬ 
tained when transitioning from the simple formulation in (13.Sp to the larger 
formulation in (EH). However, we find that, in many cases, a monolithic LP 
relaxation by itself is not able to provide tight bounds on the optimal value. 

Example 4 Let’s consider the Himmilbeau function taken from |12) . shown as 
example ID 1 in Table [TJ The POP is given by 

p{xi,X2) = {xi"^ + X2- 11)^ -b {xi +X2^ - 7)^ On [-5,5]^. (3.10) 

Solving the LP formulation (13.61) yields ps^'^'^ = —856.42 . If one used 
relaxation (EH), then we have a linear program with 324 variables and 341 
constraints, without counting the roughly 628 bounds constraints on our vari¬ 
ables. Instead, we can solve the linear program given by EH- In that case, 
we only have 25 decision variables. The matrix As will contain 200 rows and 
25 columns. Using the iterative approach, however, we just need 3 iterations 
to obtain ps^'^'^ where the matrix As contains 6 rows, in all. Thus, we achieve 
a significant reduction in the size of the LP and hence the cost of solving it. 
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However, in spite of these improvements, the objective value when us¬ 
ing (lil.8|) is —856.416. This is a very coarse lower bound on the actual optimal 
value which is p* = 0. One reason for getting a poor bound is that the consid¬ 
ered box is relatively big and that the optimal solution x* is located quite far 
from the edges. 

This motivates the Branch and Bound algorithms we are going to present in 
the next section. Before doing that, we will briefly show how one can extend 
the previous relaxations in the case of non rectangular domains. 


3.4 Extension to polyhedral and semi algebraic sets 

If /C is a bounded polyhedral set, our POP can be formulated as follows : 


minimize p{x) 
s.t X G 


a; G [0,1]” 
AqX < bo, 


(3.11) 


where Aq G and bo G K™. In fact, it suffices to compute a bounding box 

for the polyhedral set K, and then map the problem to the unit box. 

Proposition 5 Using the same notation, we build the following LP: 


PS* = minimize bs ■ zs 

s.t Os < zs < us, 


\s ■ zs = 1, 

Aszs < cs, 


(3.12) 



i<s 


Then ps* < p*, where p* is the optimal value of 


Proof The proof follows directly from the following property : 


Va; G [0,1]", jBsj{x) = x. 


i<s 


Now, If Af is a bounded semi-algebraic set, our POP can be formulated as 
follows : 


minimise p{x) 

s.c a;G[0,l]", (3.13) 


(3.13) 


gi{x) <0, V* = 1,..., TO. 


where p and gi are multivariate polynomials of degree less than 5 for all i 
1,..., TO. Then, we have the following result : 
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Proposition 6 Recall LP (13.121) below: 

PS* = minimize hs{p) ■ zg 

s.t Qs < zs < us, 

U ■zs = l, (3.14) 

Agzs < cs- 

bs{gi) • 25 < 0, Vi = 1,... ,m. 

where bs(p) = (bij{p))i<s and bs{gi) = {bi,sigi))i<s are Bernstein coefficients 
of respectively p and gi for all i = 1,... ,m. 

Then ps* < p*, where p* is the optimal value of iS.1,% . 

Proof It suffices to write polynomials gi, for all i = 1, ..., m, in the Bernstein 
basis up to the degree S and replace Bernstein polynomials using fresh variables 
zj^s for all I < S. 


4 Precision Improvements 

We will now consider three different approaches to improving our relaxation 

using the improved LP formulations proposed in this section: 

(a) We will show how further properties of Bernstein polynomials can result in 
multiaffine constraint system that can be converted back into a LP through 
dualization. However, we will see that doing so yields impractically large 
LPs. Therefore, this approach is of theoretical interest. 

(b) Next, we will consider using higher degrees 5 in our LP formulations be¬ 
yond the degrees of the original POP. However, we observe that the con¬ 
vergence is linear in |, and thus quite poor when compared to the growth 
in running times. 

(c) Finally, we will use a branch-and-bound scheme that decomposes our prob¬ 
lem domain into multiple smaller boxes, using many pruning ideas to limit 
the number of branches needed. In this context, we examine whether the 
improved LP relaxations can translate into fewer decompositions of the 
feasible region. 


4.1 Further Valid Inequalities 

We now consider techniques for adding further valid inequality constraints to 
the overall problem. As before, our goal is to ensure that the added constraints 
are affine, or can somehow be converted to an affine system of constraints. 

Adding Known Positive Polynomials: One simple approach, following 
recent developments in so-called diagonally dominant sum-of-squares is to add 
polynomials that are easy to show nonnegative such as Dij (x) : {Xi—Xj > 0 
and Eij (x) : (x^ -l-Xj)^'^ > 0, for pairs Xi, Xj to the system of constraints for de¬ 
grees d < ^ min((5i, 6j) [T]. To add such polynomials, we convert Dij and Eij to 
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the Bernstein basis, perform RLT by replacing Bernstein polynomials Bi^s{x) 
with a fresh variable Zi^s- The resulting constraints will also be added to the 
matrix {As^ cs) and possibly included in the matrix {As, cs). However, the cone 
of positive polynomials over K, is not finitely generated cone (even when we 
consider positive polynomials of bounded degrees). Therefore, an addition of 
finitely many generators cannot be useful for all problems, in general. 

Adding Multiaffine Constraints 

In this section, we briefly sketch a further approach to LP relaxations that 
involves adding multiaffine constraints and relaxing the resulting set of con¬ 
straints back to a linear program. The multiaffine constraints are given by 
product of Bernstein polynomials. Consider Bernstein polynomials pi{x) := 
Bii,Si{x),P 2 ■= ■■■ Pj '■= Bi-^Sj- 

Claim The product pi x p 2 x ■■■ x pj is of the form c(/i , ... ... ,5j) x 

_h/j.iSi-i_ sSj ; where c(/i, ... ... ,5j) is a constant coefficient given 

by the ratio of the binomial coefficients. 

This allows us to provide additional constraints in the formulation (13.611 of 
the form; 

If 1 ^/ 2 ,772 • ■ = c(/l, 5j)zi^ + ,„j^^Ki + ...+Kj 

The addition of these constraints yields a system of linear multiaffine con¬ 
straints of the following form: 

min czs 

s.t. Agzs < Cs Linear relationships between Bernstein polynomials 
zs C [Ls,Us\ bounds constraints 

Zii ''' z.i - = CiZi multiaffine equality constraints 

(4.1) 

As such, the multi-affine system above is, in fact, a nonlinear system of 
constraints. However, the following result by Ben Sassi and Girard [14], shows 
that any such system can be relaxed to yield a linear programs. 

Claim (Ben Sassi + Girard m . The multi-affine formulation in Eq. 63) 
can be relaxed to yield a linear program whose optimal value lower bounds 
that of the multi-affine system 63- 

The central idea behind Ben Sassi and Girard’s result involves writing down 
the Lagrangian L{z, p, A) involving the primal variables z and multipliers p, X 
for the equalities and inequalities in the optimization problem 63- It is noted 
that the function L is multi-afHne in z, and also that the optimal value of a 
multi-affine function in a box {Ls, Us\ is achieved at its vertices. Therefore, the 
dual is obtained as rahiz L{p,X) = min^^v L{v, p, X), where V represents the 
verticesof the box \Ls, Us\- As a result of this, the resulting LP is exponential 
in the number of variables in zs, which is already O(nl^l). 

As a result, even the addition of additional multi-afHne facts involving 
Bernstein polynomials can cause an unacceptable blowup in the problem size. 
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Table 2 Improvement of the lower bound by considering higher dimension relaxations 


&’ 

(5,4) 

(5,5) 

(6,6) 

(10,10) 

(20,19) 

(20,20) 

(^) 

-738.918 

-582.783 

-436.57 

-165.89 

-63.89 

-62.23 


4.2 Higher Degree Relaxations 

To improve the precision of the computed lower bound one can increase the 
degree of the relaxation 6. However, if we use the simpler formulations in 
Eq. (I3.5|) . then increasing S alone does not necessarily yield a better optimal 
value. 

Example 5 In the Example [H we saw that for 5 = (2,2), the optimal value 
ps^^^ = —0.5 using formulation in Eq. (0^. Now increasing the degree to 
5' = (3, 2), one can verify that = —0.59 which is a worse bound. 

However, if we used the formulation in Eq. (13.61) or the equivalent formula¬ 
tions in (13.71) and (13.81) . then it is easy to see that increasing the degree S will 
result in the addition of more constraints to the LP and thus, cannot make 
the lower bound worse. Increasing the degree of the approximation eventually 
results in tighter bounds that asymptotically converge to the globally optimal 
bound. This is motivated by the following result by Lin and Rokne [5]: 

Proposition 7 For a degree K € N", let bj K denote Bernstein coeffieients 
for a polynomial p. Then: 

Nevertheless, this convergence can be quite slow in practice. 

Example 6 Let’s consider again the POP (13.3L by increasing the degree, we 
obtain the results reported in Table[2J The results show initially large improve¬ 
ments upon increasing the degree. However, it is clear that large degrees are 
needed to approach the optimal value of p* = 0. 

This motivates us to consider the approaches developed in the previous sec¬ 
tion inside a branch-and-bound solver that recursively partitions the feasible 
region into smaller region, while lower bounding the optimal value inside each 
region using the approach considered here. In this setting, a better lower bound 
can potentially lead to fewer branches, and therefore a better performance. 


4.3 Branch-and-bound scheme 

In this section, we consider the branch-and-bound approach for solving POPs 
and integrate the improved LP formulation in Algorithm [I] into our overall 
branch-and-bound scheme. Our branch-and-bound scheme is built on top of 
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1 

2 

3 

4 

5 

6 


7 

8 


9 

10 

11 
12 


Data: Objective: p{x), constraints gi(x) > 0,. .. ,gk{^) ^ 0 domain x £ }C. 

Result: Lower bound p to the optimal value of POP. 

begin 

worklistOfBoxes := { K. } 
glbMin := +oo 

while worklistOfBoxes ^ 0 do 
B := pop (worklistOfBoxes ) 

/* Call Algorithm [J as a subroutine 

(pB,zB) := Compute Bernstein lower bound for p{x) on B 

/* Check if we have to brainch further 

/* 1. Is the value computed exact for B 

exact := Check if (pB,zB) is an exact solution to B 

/* 2. Monotonicity check. 

monotone := Check if is sign invariant over B for each Xr- 

/* 3. Heuristic termination condition, eg., box size is below 
threshold 

terminal := Check if we can terminate the branch-and-bound for B 
if monotone then 

/* 4. Create edge subproblem p,B. 
pB := branchAndBound(p, B) 
glbMin := min(glbMin, pB) 


13 

14 


else if terminal or exact then 
glbMin := min(glbMin,pB) 


15 

16 
17 


else 

/* Brauich into multiple subboxes 
splitBox {B) 

add Bi,... ^ Bk to worklistOfBoxes 


*/ 

*/ 

*/ 

*/ 

*/ 

*/ 


*/ 


Algorithm 2: Basic Branch-And-Bound scheme for solving POPs. 


previous work by Nataraj et al. m that is based on a simple formulation that 
involves finding the minimum Bernstein coefficient inside each box decomposi¬ 
tion considered by the algorithm. Additionally, their approach uses properties 
such as the vertex condition and a monotonicity condition (described below) 
to detect leaf nodes. We augment our approach directly inside their framework 
by iteratively solving LPs as described in Algorithm [TJ While solving a LP is 
more expensive than finding the minimum Bernstein coefficient, we show that 
the extra overhead is offset by our ability to consider fewer boxes. 


4-3.1 Overview of Branch and Bound Algorithm 

The main idea of the branch-and-bound (BB) algorithm is to keep subdividing 
the rectangular domain into sub-boxes until a termination condition can be 
obtained. Algorithm [5] shows the basic branch and bound scheme. It involves 
repeated decompositions of the original box /C to construct a worklistOfBoxes 
that should become empty (ideally) in order to ensure termination. 

The algorithm’s behavior and performance depends critically on three key 
operations: (a) The precise relaxation used to compure the bound for p{x) in 
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line El (b) The exactness test in lines [7] and termination check in line El and 
(c) The branching step in line 1161 

4-3.2 Exactness Test 

The exactness test is performed to infer if the current lower bound pB for p{x) 
over a given box B is in fact the optimal value. This is achieved by testing 
for the vertex condition and a monotonicity condition. The vertex condition is 
described in Corollary [2l (page [71). This is quite easy to test once we transform 
the problem from the current box B to [0,1]" using the mapping x' = T(x), 
and compute the Bernstein coefficients of p(T(x)). 

4-3.3 Monotonicity Test 

The monotonicity test (originally proposed by Nataraj et al. [H]) checks 
whether 0 ^ for r = 1,... ,n where -x.GB.li the partial derivative Pr 
w.r.t some Xr is sign invariant over B, then the global minimum of p in B 
can be obtained at one of the bounds: Xr = £r or Xr = Ur, depending on the 
sign of Pr- The derivative Pr is also expressed using Bernstein polynomials, 
where the coefficients are computed directly from the Bernstein coefficients of 
p. The monotonicity test is computed along each dimension Xr by computing 
the Bernstein coefficients of . If the polynomial is deemed monotone along 
Xr, then depending on the sign of the partial derivative, Xr is substituted by 
its lower (partial derivative is positive) or upper (partial derivative is nega¬ 
tive) bound in B. In particular, further decomposition of B is unnecessary in 
this case. However, since the global minimum may lie along a facet, we cre¬ 
ate an “edge” subproblem p by substituting Xi = li for each monotonically 
increasing variable Xi and Xi = Ui for each monotonically decreasing Xi- The 
resulting subproblem has strictly fewer variables than the original problem, 
and is solved recursively using the same branch-and-bound procedure. 

4-3-4 Termination Test 

The main termination test compares the current lower bound for the box 
B against the best upper bound p obtained by sampling feasible points in 
the original feasible region 1C. If the lower bound pB < (1 — e)p (alternatively 
pB < —e when p = 0), we do not subdivide the box further. Another approach 
to cutting off the branch-and-bound imposes a bound on the volumes of boxes 
that can be subdivided. 

4-3.5 Computing Lower Bounds 

Next, we consider the computation of lower bounds to a polynomial p{x) over 
a box B. This is a key step in our branch-and-bound scheme. We consider the 
three relaxtions defined in Eqs. dSa, (ESI) and (EH). As mentioned earlier, 
using (EH) is equivalent to computing the minimal Bernstein coefficient as 
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originally suggested by Nataraj et al. m- However, the relaxtions in (IS3]) and 
(ESI) involve solving linear programming problems that are more expensive 
when compared to hnding the smallest Bernstein coefficient. On the other 
hand, the advantage is that we obtain tighter bounds that may allow us to use 
fewer decompositions. 

As a further optimization, we build a function called “First-LP” that at¬ 
tempts to provide a lower bound for (13.511 directly without using a LP solver 
by finding a dual feasible solution for it. We rewrite (13.51) as follows: 


min b*z 
s.t. z >0 
—z > —u 
l*z = 1 


(4.2) 


wherein b is the vector of Bernstein coefficients and u represents the vector 
of upper bounds. Let us sort the Bernstein coefficents in b and without loss 
of generality we write: 

^ ^2 ^ ^ • 

Next, let bi be an index such that 6; < 0 and bi+i > 0. If &i > 0 then z* = 0 
is an optimal solution to (14.21) . On the other hand, if bN < 0, then we take 
I = N. Note that &i is the optimal value for the relaxation EU- Next, we 
choose the index q = max{* G [1, / — 1] | % — !}■ 

Lemma 3 The optimal value of (14.21) is lower bounded by max(&i,&q+i -|- 

bjUj). 

Proof We first formulate the dual to (021). Let us use the multiplier A corre¬ 
sponding to the upper bound constraints —z > u and p, corresponding to the 
equality constraint l*z = 1. The (simplified) dual LP is given as 

max /i — u* A 
s.t. pi — X < b 
A >0 

We set the dual solutions as Xj = —bj for j £ [1, g] and Xj = 0 for j > q. Finally 
we set p = bg+i. We can verify that all the dual constraints are satisfied. Thus 
our solution is dual feasible. We also note that it yields a dual objective value 
of bjUj as required. In contrast, setting p = bi and A = 0 yields 

another dual feasible solution. The rest follows by applying the standard weak 
duality theorem for linear programs. 

It is possible to provide precise conditions under which the dual feasible so¬ 
lution is in fact dual optimal, and obtain a corresponding primal optimal 
solution. The advantage of using a dual lower bound in a branch-and-bound 
scheme is that it provides an improved bound over E31) but at a reduced com¬ 
putational cost that involves sorting the Bernstein coefficeints and performing 
a linear time scan over them to identify the indices /, q which is less expensive 
than solving p.5p . For (13.6p . a lower bound is obtained by considering the 
optimal value given by First-LP, construct an associate feasible solution to it, 
then perform an iterative approach to improve this optimal value. 
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4-.3.6 Numerical Results 

The algorithms described thus far were implemented inside the MATLAB(tm) 
environment using the inbuilt linprog function for solving linear programs. 
We compare our three algorithms using a set of 18 benchmarks to evaluate 
whether the additional inequalities lead to (a) fewer boxes being examined by 
our branch-and-bound scheme and (b) overall improvement in the computation 
time. The first 16 benchmarks are collected from Nataraj et al [12] (taken in 
the same order). In addition to those, we consider two further challenging 
examples: 

— The 3-dimensional Motzkin example (1D=17) : 

p{xi,X2,X3) = Xi^X2^ +Xi'^X2'^ - 3xi^X2^X3'^ + 2 : 3 ®, R = [-0.5,0.5]^ . 

— The 4-dimensional algebraic example (1D=18): 

p { xi , X 2 , X 3 , X 4 ) = -I-2:2"^ -1-2:3'* 4-2:4"^ - 4a;ia;22;32;4 — 1, i? = [—0.1,0.1]* . 

A termination test threshold e = 10“® is fixed for computing the global 
minimum for the first 16 benchmarks. For the Motzkin example ID 17, we fix 
e = 10“^ and e = 10“^ for example ID 18 to deal with numerical issues in 
using the MATLAB’s LP solver. We expect commercial LP solvers such as 
CPLEX to provide us with more robustness. 

Table[3]shows the results obtained for the various benchmarks using the LP 
relaxations labeled 0, 1 and 2, respectively in column Ineq. These correspond 
to the LPs in (lT4l) . (IT5]) and (II3 while (1321) is computed exactly (since it 
is only given by the smallest Bernstein coefficients) whereas only lower bounds 
are computed for (1^31) and (13.61) using the results of the previous section. 
For completeness, we also report, separately, the results over the subproblems 
generated by the monotonicity tests. 

Comparing number of subdivisions: Did the use of a larger LP at each 
step yield fewer cells? From Tabled we observe that indeed the use of a larger 
LP formulation with more inequalities did lead to roughly a 10% reduction in 
the number of cells examined, especially for the larger instances. 

Comparing total time: Despite the reduction in the number of cells, 
the overall computation time for LP relaxation 2 was slightly larger. This is 
clearly due to the cost of the iterative approach (since some LPs need to be 
performed). However, for relaxation 1, the lower bound given by the First-LP 
avoid us solving LPs, which turns out to be advantageous. Indeed, the advan¬ 
tage vanishes as soon as we use an LP solver for relaxation 1, as demonstrated 
by a single example in Table |T| 

Accuracy of Results: Because of the adaptive nature of our branch-and- 
bound scheme, we obtain solutions that are consistently close to the actual 
global optima. 
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Table 3 Performance of the Cuts for benchmark problems taken from Nataraj et al m+ 
two more examples. Legend: ID: problem ID as given in \m + two more examples, Ineq.: 
the LP used for lower bounding 0: LP II3.4II . l:lower bound on II3.5II . 2:LP II3.6II . Sub. the 
number of subdivisions, Time: time taken in seconds. Cutoff: number of boxes removed 
using the cut-off test, Mono: number of boxed removed using the monotonicity test, Opt: 
the optimal value. Sub*, Cutoff*, Time*: Total number of subdivisions, number cutoff and 
time spent solving recursive subproblems. 


ID 

Ineq. 

Sub 

Time 

Cutoff 

Mono 

Sub* 

Cutoff* 

Time* 

Opt 

1 

0 

164 

1.2 

55 

62 

5 

7 

0.02 

0 

1 

1 

155 

1.1 

67 

46 

5 

7 

0.02 

0 

1 

2 

147 

2.5 

47 

61 

5 

7 

0.02 

0 

2 

0 

100 

1.1 

14 

61 

2 

6 

0.01 

-1.032 

2 

1 

97 

1.1 

15 

59 

2 

6 

0.01 

-1.032 

2 

2 

97 

3.0 

13 

61 

2 

6 

0.01 

-1.032 

3 

0 

194 

0.4 

80 

11 

3 

2 

0.00 

0 

3 

1 

176 

0.4 

85 

6 

3 

2 

0.00 

0 

3 

2 

173 

0.8 

78 

11 

3 

2 

0.00 

0 

4 

0 

1319 

4.8 

751 

4 

9 

3 

0.02 

0 

4 

1 

1199 

4.5 

684 

4 

9 

3 

0.02 

0 

4 

2 

1098 

23.2 

625 

4 

9 

3 

0.03 

0 

5 

0 

388 

3 

126 

146 

20 

22 

0.07 

-7 

5 

1 

371 

2.9 

134 

133 

17 

19 

0.06 

-7 

5 

2 

371 

3.4 

122 

145 

17 

19 

0.07 

-7 

6 

0 

784 

19.2 

371 

266 

36 

33 

0.50 

0 

6 

1 

763 

18.8 

371 

254 

34 

31 

0.47 

0 

6 

2 

707 

33.1 

318 

263 

32 

29 

0.66 

0 

7 

0 

0 

0.1 

0 

0 

0 

0 

0 

-36.713 

7 

1 

0 

0.1 

0 

0 

0 

0 

0 

-36.713 

7 

2 

0 

0.1 

0 

0 

0 

0 

0 

-36.713 

8 

0 

637 

15.6 

307 

190 

43 

40 

0.58 

0 

8 

1 

615 

15.2 

300 

181 

40 

37 

0.54 

0 

8 

2 

580 

27.3 

267 

185 

40 

37 

0.61 

0 

9 

0 

3 

0.1 

0 

2 

503 

307 

4.25 

-3.18 

9 

1 

3 

0.1 

0 

2 

498 

304 

4.25 

-3.18 

9 

2 

3 

0.1 

0 

2 

498 

304 

4.88 

-3.18 

10 

0 

1 

0.1 

0 

0 

0 

0 

0 

-20.8 

10 

1 

1 

0.1 

0 

0 

0 

0 

0 

-20.8 

10 

2 

1 

0.1 

0 

0 

0 

0 

0 

-20.8 

11 

0 

1794 

51.5 

1165 

464 

67 

73 

0.69 

-16 

11 

1 

1542 

44.4 

1050 

371 

66 

72 

0.69 

-16 

11 

2 

1525 

51.6 

989 

415 

66 

72 

0.80 

-16 

12 

0 

18 

0.1 

0 

1 

0 

0 

0.00 

-30.25 

12 

1 

18 

0.1 

0 

1 

0 

0 

0.00 

-30.25 

12 

2 

18 

0.1 

0 

1 

0 

0 

0.00 

-30.25 

13 

0 

101 

0.1 

6 

0 

0 

0 

0 

-0.25 

13 

1 

101 

0.1 

6 

0 

0 

0 

0 

-0.25 

13 

2 

101 

0.1 

6 

0 

0 

0 

0 

-0.25 

14 

0 

0 

0.1 

0 

0 

0 

0 

0 

-1.44 

14 

1 

0 

0.1 

0 

0 

0 

0 

0 

-1.44 

14 

2 

0 

0.1 

0 

0 

0 

0 

0 

-1.44 

15 

0 

118 

0.1 

7 

0 

0 

0 

0 

-0.25 

15 

1 

118 

0.1 

7 

0 

0 

0 

0 

-0.25 

15 

2 

118 

0.1 

7 

0 

0 

0 

0 

-0.25 

16 

0 

18 

0.4 

3 

0 

0 

0 

0 

-1.74 

16 

1 

18 

0.4 

3 

0 

0 

0 

0 

-1.74 

16 

2 

18 

0.4 

3 

0 

0 

0 

0 

-1.74 

17 

0 

17874 

1161 

5860 

452 

600 

404 

22.92 

0 

17 

1 

16775 

1081 

5772 

290 

600 

404 

23.23 

0 

17 

2 

16641 

1431 

5626 

452 

600 

404 

24.22 

0 

18 

0 

12684 

3636 

4080 

2422 

2616 

2240 

330 

-1 

18 

1 

12033 

3424 

4949 

1447 

2480 

2120 

314 

-1 

18 

2 

11983 

3967 

3941 

2414 

2416 

2062 

333 

-1 
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Table 4 Comparison between ’First-LP’ and Linprog performances 


ID 

Cut 

LP 

Sub 

Time 

Cutoff 

Mono 

Sub* 

Cutoff* 

Time* 

Opt 

1 

1 

‘First-LP’ 

155 

1.11 

67 

46 

5 

7 

0.02 

0 

1 

1 

Linprog 

140 

3.57 

62 

42 

5 

7 

0.07 

0 


^.5.7 Lyapunov Stability Proofs 

A standard approach to prove stability for polynomial dynamical systems is 
to find a polynomial Lyapunov certificate which consists on a positive definite 
function decreasing along the trajectories inside a region of interest. More pre¬ 
cisely, let fo be a polynomial candidate Lyapunov function, V its derivative 
and R the region of interest taken as a rectangle containing zero (the equilib¬ 
rium point). To verify the asymptotic stability of the equilibrium, we should 
verify that: minfo(x) > 0 and min—t7(a;) > 0. 

The advantage while solving POPs arising from Lyapunov function synthe¬ 
sis problems is that a global minimum is known in advance. In fact since usually 
y(0) = P(0) = 0, then we already know that zero is the global minimum of a 
true Lyapunov function. Therefore, a good branch-and-bound decomposition 
scheme for this problem decomposes around the equilibrium to maximize the 
opportunity for exact relaxations m- 

To show the efficiency of the zero decomposition, we consider 9 Benchmarks 
given in our earlier work m, taken in order. The goal is to verify that the 
candidate functions are indeed Lyapunov functions. In all these examples, the 
region of interest is R = [—1,1]". We propose to check the validity of these 
results by computing pv* and Py* which are lower bounds on V and —V inside 
R using the smallest Bernstein coefficient (pv*(0), py*(0)) and relaxation (13.51) 
(py*(l),Py*(!)). We report in Table [5] the results we obtained where stability 
is said verified once a precision of I0“® is reached. In the appendix we give 
a detailed description of the Benchmarks, the Lyapunov function and their 
associated Lie derivatives. 


Table 5 Proving bounds on Lyapunov functions and their derivatives. Legend: EX - ID 
of the example taken from Ben Sassi et al. |15| . PyU)- Lower bounds to optimal value 
obtained by using LP relaxation id j, Pydot^^y' Lower bounds on optimal value of Lyapunov 
derivative. 


EX 

Pv*(0) 

Pv*(l) 

Pvdoti^) 

PvdntW 

Stability 

1 

-9.2 X 10"^^ 

0 

-5.4 X 10“^^ 

-3.5 X lO"^'’ 

✓ 

2 

-1 

-0.0625 

-6.2 X 10-1^ 

0 

X 

3 

-5.4 X lO"^*^ 

0 

-2.9 X lO"^*^ 

0 

✓ 

4 

-2.7 X 10“^ 

0 

-8.7 X lO-^*^ 

-1 X 10"^“ 

✓ 

5 

-1.3 X lO"^*^ 

0 

-3.7 X 10"^^ 

-1.4 X lO”^"* 

✓ 

6 

-3.4 X 10-^^ 

-3.5 X lO-^’^ 

-6.9 X lO-^^’' 

-4.1 X 10-^^” 

✓ 

7 

-1.5 X lO"^*^ 

-2.1 X lO-^'’' 

-9.5 X 10“^^ 

-3.8 X 10“^^ 

✓ 

8 

-10.9788 

-10.9788 

-7.9 X 10“^ 

-3.3 X 10-y 

X 

9 

-9.9 X lO"^*^ 

-1.7 X lO-^’" 

-3.8 X 10“^^ 

-3.3 X 10-^^ 

✓ 
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5 Conclusions 

We present a novel approach to deal with polynomial optimization problems 
(POPs) by relaxing them to bigger size linear programs. The key idea is to use 
Bernstein polynomials in order to build LPs that can handle many of the re¬ 
lations between non linear terms missed because of the linearization process. 
Contrarily to the standard RLT approach, the given LPs are easily imple- 
mentable since only a Bernstein framework is needed (coefficients, bounds and 
change of variable). Thanks to the properties of Bernstein polynomials , tighter 
bounds than RLT are obtained and various techniques to improve the preci¬ 
sion of these bounds are given. We show that our relaxations can be used to 
improve the Branch and Bound scheme given by Nataraj |12j . The main draw¬ 
back faced in the latter case was the extra cost of solving LPs. We already find 
a way to avoid this for our first linear relaxation but not for the more precise 
one. This is definitely a first goal future work. Also, we manage to extend our 
Brand and Bound algorithms in the case of semi algebraic constraints. 
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6 Appendix 

Benchmark #1: 

dx 

dt 

dt 


Consider the two variable polynomial ODE: 
= -12.5x + 2.5x^ + 2.5y2 + IQx^y + 2.5j/^ 

= -y- y^- 


Lyapunov function : 

2x^2 5y'2 

Lyapunov derivative function : 

40x"3y+10x"3 —50x"2+10xy "3+lOxy'2 —lOy "3 —lOy "2 

Benchmark ^^2: Consider the two variable polynomial ODE: 

^ = 6.933333a;^ + 4.566667a;^ - 21.5x. 
at 

^ = 6.933333a;^ + 0.4a;^y + 2.066667x^ + xy^ + 0.6x1/ — 9x — y"^ — y. 
dt 


Lyapunov function : 

5x"2—4xy+5y" 2 

Lyapunov derivative function : 

41.6x"4+40x^3y + 37.4000x^3 -179x‘2+10xy"3+10xy‘2 -10y^3 -10y^2. 


Benchmark #3: 

dx 

dt 

dt 


Consider the two variable polynomial ODE: 
= —1.5x — x'^ + 0.5x1/ + O.by^ — 2x^ + x^y. 

= -0.5y. 
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Lyapunov function : 

5x^2+5y^2. 

Lyapunov derivative function : 

—20x"4+10x"3y—10x*3+5x"2y—15x*2+5xy"2—5y " 2. 

Benchmark #4: Consider the two variable polynomial ODE: 

dx Q 

— = —2x^ — 0.5xy — 0.5a;. 
dt ^ 

^ = 0.25xy^ - 0.125xy + 0.25^^ - 0.4125y. 
dt 


Lyapunov function : 

5x^2+5y "2. 

Lyapunov derivative function : 

-20x^4 -5x^2y -5x'2 + 2.5xy'3 - 1.25xy ^2 + 2.5y'3-4.125y ^ 2 . 
Benchmark #5: Consider the three variable polynomial ODE: 
dx 


dt 

dt 

dz 

dt 


= —2x^ — 0.5xy — 0.5a; — — z"^. 

= 0.25a;2/2 - Q.l2?>xy + 0.25^^ - 0.4125|/. 


= —z — z. 


Lyapunov function : 

5x^2+5y'2+5z '2. 

Lyapunov derivative function : 

-20x'4-5x'2y-5x^2 + 2.5xy^3-1.25xy'2-10xz'3-10xz^2 + 2.5y^3-4.125y'2-10z^3 

Benchmark #6: Consider the three variable polynomial ODE: 

dx 
dt 

dt 
dz 
dt 


= -0.5x^y + 0.5a;3z2 - 5x^ + y^ - y"^ + yz'^ - z^. 
= 0.25y2 - 0.25y. 

= yz"^ + z'^ - 2z^. 


-lOz "2. 
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Lyapunov function : 

1.9150x^4+5x~3+5x' 2y'2+5x^2z'2 + 3.9396x^2-2.5409xy^ 3 + 2.5409xy^ 2 
+5y"4+5y'3+5y^2z'2+5y"2+5zM+5z^3+5z'2 . 

Lyapunov derivative function : 

-3.8300x‘6y + 3.8300x‘6z^2-22.98x^6-7.5x^5y + 7.5x'5z^2-45x'5-5x^4y^3+5x'4y^2z'2 
-30x^4y‘2-5x^4yz^2-3.9396x^4y+5x‘4z'4-26.0604x^4z‘2-23.6375x^4 + 7.66x'3y^5 
-6.3895x‘3y^4-1.2705x^3y'3z^2 + 6.3524x'3y^3+1.2705x^3y^2z^2-7.6228x'3y'2 
+ 7.66x~3yz"4-7.66x'3z^4+15x^2y"5-15x'2y'4 + 2.5x'2y'3-2.5x‘2y^2+10x‘2yz'5 
+15x"2yz *4+10x"2z "5 —35x"2z *4+10xy "7 —lOxy '6+lOxy"5z * 2 + 7.8792xy "5 —lOxy"4z "2 
— 9.7849xy'4+lOxy"3z "4 + 3.1762xy "3 —lOxy"2z '4 — 1.2705xy "2+lOxyz "6 + 7.8792xyz "4 
-10xz^6-7.8792xz^4-2.5409y'8 + 5.0819y'7-2.5409y"6+5y'5-2.5409y"4z'4-1.25y"4 
+ 10y^3z^5 + 5.0819y'3z^4 + 2.5y^3z^2-1.25y'3+10y^2z"5-22.5409y^2z^4-2.5y"2z^2 
-2.5y'2+20yz'7+15yz'6+10yz'5+20z'7-25z‘6-20z'5-20z '4. 


Benchmark #7: Consider the three variable polynomial ODE: 


— = -O.Bx^p + 0.5a;^2;^ - + y'^z + y^ - yz^ + yz"^ + z^ - z'^ 

at 

^ = 0.5y^z — 0.5?/^ — 2y 

dz 2 2 

— = -yz +yz + z - z 


Lyapunov function : 

-1.2500x^4+1.6667x"3+5x'2y'2+5x^2z"2+5x'2+5xy^3+5xy^2z+5xy"2 + 1.0921xz^3- 
1.0921 xz"2+5y*4+5y "3z+5y*3+5y *2z"2+5y"2z+5y "2 + 0.4638yz~3—0.4638yz *2 
+5z"4+5z"3+5z "2. 

Lyapunov derivative function : 

2.5x"6y —2.5x"6z"2+5x"6 —2.5x"5y + 2.5x"5z "2 —5x"5—5x"4y "3+5x "4y "2z"2 —10x"4y"2 
—5x"4yz "2—5x"4y+5x "4z"4 —5x"4z"2 —10x"4—5x"3y"4z—7.5x"3y"4+2.5x"3y"3z"2—2.5x"3y"3z 
— 7.5x"3y"3 + 2.5x"3y"2z "3 + 2.5x"3y "2z"2—5x"3y"2z—5x"3y "2 + 4.454x"3yz "3—4.454x"3yz "2 
+ 0.546x"3z"5-0.546x"3z"4-6.0921x"3z"3 + 6.0921x"3z"2+5x"2y"4z+5x"2y"4+5x"2y"3z 
—5x" 2y "3 —20x" 2y "2 —15x "2yz"3+15x"2yz"2 + 15x"2z"3 —15x"2z "2 + lOxy" 6z+10xy" 6 
+10xy "4z "3+lOxy "4z " 2 + 17.5xy "4z +2.5xy "4 —lOxy"3z "3+lOxy"3z "2+5xy"3z—35xy"3 
+ 10xy "2z"3—5xy"2z"2 —25xy"2z—20xy"2 —lOxyz "5 + 6.7237xyz "4—4.5395xyz "3 
+ 7.8158xyz"2 + 10xz"5-6.7237xz"4 + 4.5395xz"3-7.8158xz"2+5y"7z+5y"7+5y"6z"2 
+10y"6z+5y"6+10y"5z —lOy "5 + 1.0921y"4z"4—5y"4z "3 + 6.4079y"4z"2+5y"4z 
-47.5y"4-5y"3z"4+10y"3z"2-30.0000y"3z-35y"3 + 3.8406y"2z"4 + 11.8551y"2z"3 
-30.6957y"2z"2-25y"2z-20y"2-1.0921yz"6-17.8158yz"5 + 5.2992yz"4 + 1.7535yz"3 
+ 11.8551yz"2 + 1.0921z"6 + 17.8158z"5-3.9079z"4-5z"3-10z "2. 
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Benchmark 9 ^ 8 : Consider the three variable polynomial ODE: 

^ = -0.5x^y + 0.5x^z‘^ - + y'^z + - yz^ + iyz"^ + 2 ^ - 3z^ 

at 

§ = /^-/_2y3_,3^3^2 

at 


Lyapunov function : 

2.3519x^5-1.0449x'4y + 0.3429x‘4z+5x"4 + 4.4496x^3y‘2-2.5x‘3y+5x^3z'2+5x'3+5x^2y'3+5x^2y^2 
+5x"2yz"2+5x*2z"3+5x"2z"2+5x"2+3.6461xy"4+5xy"3+5xy"2z"2+5xyz"3+5xyz"2+5xz"4+5xz"3 
+ 2.5863xz^2-0.4325y'5 + 4.9487y^4z+5y'4+5y^3z^2-4.0594y^3+5y‘2z'3+5y'2z'2+5y^2+5yz'4 
+5yz'3 + 4.5809yz'2 + 1.8568z^5+5z^4-0.7475z'3+5z ^2. 

Lyapunov derivative function : 

-5.8798x'7y + 5.8798x*7z^2 -11.7597x"7 + 2.0899x"6y ^2-2.0899x"6yz'2-0.6857x^6yz-5.8203x^6y 
+ 0.6857x^6z‘3+10x^6z^2 -1.3715x^6z-20x‘6-6.6743x'5y ^3 + 6.6743x^5y ^2z ^2-9.5986x'5y “ 2 
-11.25x"5yz^2 + 7.5x^5z'4-7.5x'5z‘2-15x^5 + 10.7148x'4y^4z + 7.8046x^4y'4+5x^4y^3z‘2 
-12.9101x^4y'3-10x^4y'2+5x'4yz'4-16.7597x'4yz'3 + 20.279x"4yz'2 

-5x‘4y+5x^4z^5+5x^4z'4 + 2.8046x^4z^3-43.071x^4z^2-1.0286x'4z-10x'4 + 4.7194x^3y'5z 
-14.9019x^3y'5 + 3.1948x'3y'4z^2 + 18.8712x^3yMz-1.4443x^3y'4+2.5x^3y'2z'4 
+ 1.6797x"3y^2z"3-20.0392x^3y"2z^2 + 2.5x'3yz'5-1.3715x^3yz"4-36.4644x'3yz^3 
+ 92.9436x'3yz'2 + 2.5x'3z"6 + 2.5x“3z'5-2.3357x^3z"4 + 23.3865x'3z‘3-100.0863x^3z^2 
+ 28.3487x'2y^6z-1.6513x^2y'6 + 2.5008x^2y^5z-47.5x'2y^5+20x~2y^4z"3+10x‘2y^4z^2 
+ 14.9998x'2y'4z-5x'2y"4-13.3487x^2y'3z^3 +30.046x^2y'3z'2 + 5.8514x"2y'2z ^3 

-17.5460x'2y'2z'2-15x'2yz'5+45x'2yz'4-22.5024x'2yz'3 + 67.5x'2yz'2+10x'2z'5 
-15x'2z'4-20x'2z'3-75x'2z“2 + 24.5844xy'7z-4.5844xy'7+25xy'6z-34.1687xy'6+20xy'5z'3 
—30xy'5+15xy'4z'4+lOxy'4z'3+15xy'4z'2+lOxy'4z+10xy'4 — 24.5840xy'3z'3 + 33.7529xy'3z'2 
— lOxy'2z'5+30xy'2z'4+5xy'2z'3 —15xy'2z'2 —lOxyz'6+20xyz'5+45xyz'4—45xyz'3+5xz'6+10xz'5 
-60xz'4-29.8273xz'3-45.5180xz'2 + 1.4833y'8z+5.8088y'8+19.7947y'7z'2 + 5.205y'7z 
-10.6745y'7+20y'6z'3-10y'6z'2-51.7682y'6z-27.8218y'6+15y'5z'4 + 6.3539y'5z'3 
-24.0617y'5z'2 + 10y'5z + 14.3566y'5 + 10y'4z'5+10y'4z'4-12.0244y'4z'3-19.4723y'4z'2 
-14.8461y'4z-20y'4-5y'3z'5-14.795y'3z'4 + 44.3852y'3z'3 + 5.8381y'3z'2-5y'2z'6 
+60y'2z'4-22.8216y'2z'3-66.5347y'2z'2-5yz'7+5yz'6 + 42.4137yz'5-22.2410yz'4 
-45.8383yz'3 + 2.5148yz'2 + 9.2838z'6-9.8460z'5-56.2589z'4 + 16.7274z'3-30z '2. 

Benchmark #9: Consider the three variable polynomial ODE: 

^ = 0.05x^yz + O.Obx^y - t)mx^z - 0.05a;2 + Q.Qbxyz + O.OSxy - 0.05x2 - 0.05x + 0.125^^2 - 0.125^^ 

+ 0.125y2z - 0.125J/2 + 0.27/2® + 0 . 2 y 2 ‘‘ - 0 . 22 ® - 0 . 22 ^ 

^ = 0.1257/^2 - 0.1257/2 + 0.1257/2 - 0.125// + 0.22® + 0.22^ 
at 

^ = -o.u=-o.u 

dt 
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Lyapunov function : 

2.5x'2 + 2.5y^2+5z"2 

Lyapunov derivative function : 

0.25x"3yz+0.25x"3y —0.25x*3z—0.25x*3 + 0.25x"2yz+0.25x"2y —0.25x"2z—0.25x"2 + 0.625xy"3z 

— 0.6250xy*3 + 0.6250xy"2z — 0.6250xy"2+xyz"5+xyz"4—xz"5—x+z "4 + 0.625y"3z—0.625y"3 + 0.625y"2z 

— 0.625y"2+yz"5+yz"4—z"3—z " 2 . 



